%pylab inline
import matplotlib.pyplot as plt
#import collections
Jenő Sólyom: Fundamentals of the Physics of Solids, Volume 1: Structure and Dynamics (A modern szilárdtestfizika alpajai I. A szilárd testek szerkezete és dinamikája)
See: page 350. Eq. (11.2.56) for simple cubic structure including 1st and 2nd neighbor interactions
# Az abra kimentesehez az alabbiakat a plt.show() ele kell tenni!!!
#savefig('fig_rainbow_p2_1ray.pdf'); # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps'); # Abra kimentese
# Abra es fontmeretek
xfig_meret= 12 # 12 nagy abrahoz
yfig_meret= 7 # 12 nagy abrahoz
xyticks_meret= 15 # 20 nagy abrahoz
xylabel_meret= 20 # 30 nagy abrahoz
legend_meret= 21 # 30 nagy abrahoz
def Dm_op(qv,K1,K2,delta1,delta2):
## q-points are in units of 2 pi/a, qv
Dm1= array(zeros((3,3)))
for d in delta1:
Dm1 = Dm1 + outer(d,d)/dot(d,d)*(1-exp(1j*2*pi*dot(qv,d)))
Dm2= array(zeros((3,3)))
for d in delta2:
Dm2 = Dm2 + outer(d,d)/dot(d,d)*(1-exp(1j*2*pi*dot(qv,d)))
Dm= K1*Dm1 + K2*Dm2
return Dm
# simple cubic
qv = [1.2,1.2,0]
K1=2
K2=1
deltaSC1= array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]),
array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])
deltaSC2=array([array([0, 1, 1]), array([0, 1, -1]), array([0, -1, 1]),
array([0, -1, -1]), array([1, 0, 1]), array([1, 0, -1]),
array([1, 1, 0]), array([1, -1, 0]), array([-1, 0, 1]),
array([-1, 0, -1]), array([-1, 1, 0]), array([-1, -1, 0])])
real(Dm_op(qv,K1,K2,deltaSC1,deltaSC2));
see page 354, Fig. 11.10.
K1=1
K2=0.5
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 3.0;
## k-points are in units of 2 pi/a
## high symmetry points in the Brillouin zone
## FONTOS: a fenti abra szerinti pontokat kell venni,
## a szomszedos spec. k pontokat kell venni.
Gp = array([0,0,0])
Xp = 0.5*array([0,1,0])
Rp = 0.5*array([1,1,1])
Mp = 0.5*array([1,1,0])
## G-X-M-G-R line
specpoints = array([Gp,Xp,Mp,Gp,Rp])
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Mp):"$M$",
tuple(Rp):"$R$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for kv in klist:
enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaSC1,deltaSC2))))
enlist= array(enlist)
eigszam= len(enlist[0,:])
specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,eigszam):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "phonon for sc, $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)
grid();
#savefig('Fig_phonon_simple_cubic_1-2_neighbors.pdf'); # Abra kimentese
## fcc lattice:
K1=5
K2=0.5
Npoints = 500; ## hany pontbol keszul a gorbe
ymax = 3.5;
## k-points are in units of 2 pi/a
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])
## G-X-W-L-G-K-U-X line
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Up,Xp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp,Wp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp]) # Solyom konyvben: II. kotet, 20.1. abra
# Az U,K pontok egybe!!!!!
specpoints = array([Gp,Xp,Wp,Xp,Up,Kp,Gp,Lp])
## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!
### kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
if not i==4: ########## Az U,K pontok egybe!!!!!
# a i== 4 jeloli, hogy hol van az Up pont,
# azaz az Up indexe a 'specpoints' listaban.
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
deltaFCC1= array([array([0, 1/2, 1/2]), array([0, 1/2, -1/2]), array([0, -1/2, 1/2]),
array([0, -1/2, -1/2]),
array([1/2, 0, 1/2]), array([1/2, 0, -1/2]), array([-1/2, 0, 1/2]),
array([-1/2, 0, -1/2]),
array([1/2, 1/2, 0]), array([1/2, -1/2, 0]), array([-1/2, 1/2, 0]),
array([-1/2, -1/2, 0])])
deltaFCC2=array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]),
array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])
enlist = []
for kv in klist:
enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaFCC1,deltaFCC2))))
enlist= array(enlist)
eigszam= len(enlist[0,:])
specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,eigszam):
plot(kk,enlist[:,ii],color='r')
specNev = ["$\Gamma$","$X$","$W$","$X$","$K$","$\Gamma$","$L$" ]
#specpoints = array([Gp,Xp,Wp,Xp,Up,Kp,Gp,Lp])
i = 0
for xc in specpos:
axvline(x=xc,color='b')
#text(xc,-0.5,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
text(xc,-0.5,specNev[i],fontsize=xylabel_meret)
i=i+1
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
#ylim(0,ymax)
cim = "phonon for fcc (Au) $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)
grid();
#savefig('Fig_phonon_fcc_1-2_neighbors.pdf'); # Abra kimentese
# bcc lattice
K1=1
K2=0.5
deltaBCC1 = array([array([1/2, 1/2, 1/2]), array([1/2, 1/2, -1/2]),
array([1/2, -1/2, 1/2]), array([1/2, -1/2, -1/2]),
array([-1/2, 1/2, 1/2]), array([-1/2, 1/2, -1/2]),
array([-1/2, -1/2, 1/2]), array([-1/2, -1/2, -1/2])])
deltaBCC2 = array([array([1, 0, 0]), array([-1, 0, 0]), array([0, 1, 0]),
array([0, -1, 0]), array([0, 0, 1]), array([0, 0, -1])])
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 3.0;
## k-points are in units of 2 pi/a
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Hp = 0.5*array([0,1,0])
Pp = 0.25*array([1,1,1])
Np = 0.25*array([1,1,0])
## G-H-N-G-P-H line
specpoints = array([Gp,Hp,Np,Gp,Pp,Hp])
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Hp):"$H$",tuple(Pp):"$P$",
tuple(Np):"$N$"}
figsize(xfig_meret,yfig_meret)
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for kv in klist:
enlist.append(sqrt(eigvalsh(Dm_op(kv,K1,K2,deltaBCC1,deltaBCC2))))
enlist= array(enlist)
eigszam= len(enlist[0,:])
#specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,eigszam):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$\omega(q)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
#ylim(0,ymax)
cim = "phonon for bcc, $K_1$ ="+str(K1)+ ", $K_2$ ="+str(K2)
title(cim,fontsize=xylabel_meret)
grid();
#savefig('Fig_phonon_bcc_1-2_neighbors.pdf'); # Abra kimentese